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We describe a novel switching algorithm based on a "reverse" Monte Carlo method, in which the 
potential is stochastically modified before the system configuration is moved. This new algorithm 
facilitates a generalized formulation of cluster-type Monte Carlo methods, and the generalization 
makes it possible to derive cluster algorithms for systems with both discrete and continuous de- 
grees of freedom. The roughening transition in the sine-Gordon model has been studied with this 
method, and high-accuracy simulations for system sizes up to 1024 2 were carried out to examine the 
logarithmic divergence of the surface roughness above the transition temperature, revealing clear 
evidence for universal scaling of the Kosterlitz-Thouless type. 

PACS numbers: 05.10.Ln, 05.50.+q, 64.60.Ht, 75.40.Mg 



Large-scale Monte Carlo (MC) simulations are often 
plagued by slow sampling problems. These problems are 
especially severe in systems near the critical point or in 
those with strong correlations. Slow sampling problems 
manifest themselves as poor scaling of the dynamic relax- 
ation time with the system size, making large-size sim- 
ulations extremely slow to converge. The cause of these 
problems is that most MC simulations are based on lo- 
cal moves, and when the correlation length of the system 
grows or as relaxation modes of the system become heav- 
ily entangled, local moves become increasingly inefficient. 
But if nonlocal MC moves are used their acceptance 
ratios are often found to be exceedingly low when system 
correlations are strong. 

One way to circumvent these problems was suggested 
by Swendsen and Wang Q , who devised a clever scheme 
where large-scale nonlocal MC moves may be constructed 
to achieve high sampling efficiencies by exploiting certain 
geometric symmetries in the system. This algorithm led 
to a marked reduction in the dynamical scaling exponent 
in the 2-dimensional Ising model near criticality. Since 
the nonlocal moves in this algorithm update a large num- 
ber of degrees of freedom at the same time, the Swendsen- 
Wang method and others inspired by it are also often 
referred to as "cluster Monte Carlo" methods. 

Since Swendsen and Wang's paper in 1987, many 
cluster-type MC alg orithms have appeared SJj, H H 

0, H S El EH, El Ei El El El El EaliJ . But 

the success of cluster MC has not been universal because 
the proper cluster moves needed seem to be highly de- 
pendent on the system, and efficient cluster MC meth- 
ods have been found for only a small number of models 
[I S, H 0, EH, E3, El, ES El so far. The difficulty of 
formulating a generalized MC method that could work 
for any system seems to be associated with the appar- 
ent geometric nature of the cluster-type MC methods - 
all existing cluster MC methods have been derived in 
one way or another by using certain geometric features 
of the system. For example, in the original Swendsen- 



Wang formulation a mapping between the Ising model 
and the percolation model originally described by For- 
tuin and Kasteleyn [2(| was exploited to effect cluster 
spin flips. In other models, the requisite mapping is not 
always obvious, making cluster MC methods difficult to 
implement for general systems. 

In this letter, we will show that the derivations of clus- 
ter MC methods do not have to be based on geometric 
features of the systems. Instead, they may be more con- 
veniently formulated based on algebraic features of the 
system potential V(C). We will exploit this algebraic for- 
mulation and suggest a way to generalize cluster Monte 
Carlo methods to systems with any potential. 

We focus on classical systems with partition function 
Z = Tr e~ v ^ c >, where V(C) is the potential in units of 
the temperature T. Acceptable Monte Carlo methods 
to sample the system configurations can be constructed 
using any transition probabilities W(C — > C) as long as 
the detailed balance condition 



(1) 



is satisfied. Conventional MC methods such as Metropo- 
lis [2l| accomplishes this in two steps: a trial move is 
made from C to C with some transition probability, and 
the move is then accepted or rejected according to an 
acceptance probability based on V(C), V(C) or both, so 
that the composite process satisfies Eqn.(l). This way 
of constructing the Markov chain - trial moves followed 
by acceptance/rejection - has been the accepted "stan- 
dard" method for doing MC since the inception of the MC 
method j2^|. Other MC methods do exist, such as the 
heat bath algorithm [23| , which follow alternative strate- 
gies, but by far the standard method is conceptually the 
simplest and most convenient in practice. 

In the Monte Carlo method we are proposing, we will 
reverse the order of the two steps in the standard method. 
That is, we will first determine an acceptable way to mod- 
ify the potential V and then find a transition C — > C that 
is consistent with the new potential. To our knowledge, 
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the basic elements of this "reverse MC" idea were first 
suggested by Kandel et al. [4], who used it to stochasti- 
cally remove interaction terms from the system's poten- 
tial in an Ising model to arrive at an alternative deriva- 
tion of the Swendsen-Wang method. The formulation of 
Kandel et al. was limited to discrete systems like the 
Ising model. In the following, we will show how the re- 
verse MC idea may be formulated more generally for any 
system, discrete or continuous, and how it may then be 
used as a framework to construct generalized cluster al- 
gorithms. 

Consider a system with potential V = Uj + V r , con- 
sisting of a number of "interaction terms" Vi plus a "resid- 
ual" V r . These interactions may be bonds between parti- 
cles, interactions of the particles with a field, or any other 
additive terms in V. We consider replacing each interac- 
tion term Vi by some 'pre- selected Vi with a "switching" 
probability 



St(C) 



where Avj — Vj — Vi 
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and Av* = maxc Vi(C). 



The outcome of the switches defines two complemen- 
tary sets of interactions - the switched ones a and the 
unswitched ones a. Using the outcome of the switches, 
we define a stochastically modified potential V as follows: 



(3) 
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with Vi — vi — ln(l — Si). An MC pass starts with an 
attempt to switch every interaction Vi to the new x>i us- 
ing the Si defined in Eqn.(2). If the switch is successful, 
the interaction is replaced by Vi. If not, the interaction 
is replaced by another interaction Cj. This is followed 
by an update in the configuration of the entire system 
from C to C using a transition probability W(C — > C) 
that satisfies detailed balance on the modified potential 
V. This constitutes one pass. The move from C to C can 
of course be carried out using any conventional MC move 
that satisfies detailed balance on the modified potential. 
But the reverse formulation of the MC method now offers 
possibilities that were previously unavailable to conven- 
tional MC methods — if a simple scheme can be devised 
to update the configuration of the entire system on the 
stochastically modified potential, one can envision de- 
signing global moves for the system to accelerate its sam- 
pling, and our freedom in choosing the i>i can be actively 
exploited to facilitate this. Within this context, the origi- 
nal formulation of Kandel et al. corresponds to switching 
Vi to Vi = 0, i.e. simplifying the potential by deleting in- 
teractions from it. Kandel et al. showed that for the Ising 
model they could easily construct global moves on this 
stochastically simplified potential and their formulation 
regenerates the Swendsen-Wang method. But compared 
to the deletion formulation of Kandel et al, the switch- 
ing implementation of the reverse MC method now offers 



a much wider set of possibilities because the form of the 
"switch to" interactions is completely arbitrary. Whereas 
previously there may not be an obvious way to globally 
update the configuration of the system on the original po- 
tential, with the proper choices for Vi large-scale moves 
may now become possible on the stochastically modified 
potential. Indeed, we have shown that the switching idea 
may be used to formulate a cluster MC algorithm for a 
Lennard- Jones fluid [2~i ]. 

Equations (2), (3) and the transition probability W 
define the switching algorithm. To prove detailed bal- 
ance Eqn.(l) for the switching algorithm, it is suffi- 
cient to treat a case where there are only two interac- 
tion terms. Extension to any number of interactions 
is straightforward. Starting with C, with two interac- 
tion terms v\ and i>2, there are four possible outcomes 
from the switch: I. both 1 and 2 are switched, which oc- 
curs with probability Pi — Si(C)S 2 {C), II. 1 is switched 
and 2 is unswitched, with P u = Si(C)[l - S 2 (C)], III. 
1 is unswitched and 2 is switched, with Pm = [1 — 
Si(C)]S2(C), and IV. both 1 and 2 are unswitched, with 
Prv = [1 - Si(C)} [1 - S 2 (C)]. After the switch, an up- 
date C — > C is made with a transition probability W 
that satisfies detailed balance on the modified potential 
V defined in Eqn.(3). Each of the four channels will 
have a different W: W h Wn, etc., and W(C -> C) in 
Eqn.(l) is the sum PiWi + ffrWii + fin Win + PiyWrv 
over all four channels. For the reverse transition, we 
start with C and consider switching Ui(C') — > vi(C) and 
v 2 (C) — > v 2 (C'). Again there are four possible outcomes 
and we call these scenarios I', II', III' and IV' as for 
the forward transition. W(C — ► C) in Eqn.(l) is again 
the sum P v Wy + PwWiv + fin' Win' + Prv'Wiv Using 
the choice of S and V in Eqs.(2) and (3), it is easy to 
show that detailed balance is obeyed along each chan- 
nel, i.e. PiWi = PyWy, PxiWu = PiyWw, etc. Of 
course, detailed balance only requires the total W to sat- 
isfy Eqn.(l), and it is possible to choose alternate forms 
of S and V to do that, which may provide further flexi- 
bilities. 

In the rest of this letter, we will illustrate the effec- 
tiveness of the switching implementation of the reverse 
MC method, and show how it can be used to easily de- 
rive a cluster MC method in a system with continuous 
degrees of freedom. Previously, it has been extremely dif- 
ficult to design cluster MC algorithms for systems with 
continuous degrees of freedom. The few that have been 
reported to date 0, [|| 0, E , 14 , 19[ were mainly based 
on embedding discrete degrees of freedom into continu- 
ous ones. The only exception is the recent discovery of a 
geometric MC algorithm by Liu and Luitjen [l^ where 
they formulated a rejection-free MC method to sample 
the Lennard- Jones fluid at its critical point. 

The switching algorithm we have proposed makes the 
process of deriving cluster-type MC methods much more 
straightforward compared to those based on geometric 



features of the system. We will illustrate this using 
the sine-Gordon model, which can be used to study the 
roughening transition on 2-dimensional surfaces. The 
sine-Gordon (SG) model has the potential 



^sg = T- 
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where <f>i are continuous variables on a 2-dimcnsional 
square lattice, the second sum is over all sites and the 
first sum is over all nearest-neighbor pairs. The SG 
model is often considered to be a coarse-grained ver- 
sion of the discrete Gaussian (DG) model with poten- 
tial Vdg = j) \hi ~ hj\ 2 , where hi are integers. 
The DG model can in turn be mapped directly onto 
the Coulomb gas model [2a ]. and as a result, the SG 
model should belong in the same universality class as the 
Kosterlitz-Thouless (KT) transition [H,[27|]. 

Roughening is expected to be a weak transition. The 
only easily discernible divergence is exhibited in a loga- 
rithmic dependence of the surface roughness a 2 = (\<j>i — 
(4>)\ 2 ) on the system size L at the roughening tempera- 
ture Tr. Below Tr, a 2 is expected to approach a finite 
value as L — > oo. In addition to this, since the divergence 
is slow, large lattice sizes are needed to reach the scaling 
limit. All of these features of the SG model make it hard 
to accurately study the roughening transition using MC 
simulations. Previous simulations have been limited to 
small systems [13, 0, M, E M, El . 

In order to locate Tr and study the scaling behavior 
at the roughening transition, we make use of the switch- 
ing algorithm of the reverse MC method proposed above. 
The essential difficulty in treating the SG model is due 
to the nonlinear cosine terms in the potential in Eqn. (4) . 
If these nonlinearities could be removed, the residual po- 
tential becomes a simple Gaussian and we could move 
the system configuration effectively using uncoupled sur- 
face modes. With this in mind, we separate the potential 
into two parts and treat the cosine terms as interactions 
Vi = — T cos 4>i and the harmonic part as the resid- 
ual V r . Each of the interactions is switched to a uni- 
form potential v t = -T~ x with Si = el 1 "™ 8 ^]/ 7 . After 
the switches, a number of <pi would have effectively lost 
their couplings to the cosine potential, while the rest have 
their interactions with the cosine potential replaced by 
Vi = — ln[e COS( ^/ T — e~ 1 l T \. In the ensuring MC move, 
we can update the unswitched 0, which are now cou- 
pled to the replacement interactions using conventional 
methods, but try to formulate an update scheme where 
the rest of the <pi, now forming a constrained Gaussian 
field, may be updated globally. A Gaussian field sub- 
ject to linear constraints is still Gaussian, and in prin- 
ciple we can diagonalize the potential to obtain all the 
normal modes and then move each one independently. 
This problem is the subject of fracton dynamics and has 
been studied previously [3^|. However, the cost of ob- 
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FIG. 1: (a) Surface thickness a as a function of the log 
of lattice size L for different temperatures T. (b) Expanded 
view of (a) for several temperatures near Tr shifted vertically 
to coincide at L = 64. Dashed line is the expected KT slope 
at Tr, showing that Tr is slightly above T = 25 but below 
26. 



taining all the normal modes of the constrained surface 
and their frequencies will grow rapidly with the size of 
the lattice and will only be feasible for small-size simula- 
tions. Since the scaling limit in the SG model can only 
be reached with large system sizes, we will need an al- 
ternative method. The method we have used to update 
the constrained Gaussian fields is based on the method 
of Hoffman and Ribak 



34(. Since the statistics of the 



fluctuations of a Gaussian field from its mean is indepen- 
dent of the value of the mean field, the fluctuations from 
a free Gaussian field can be transferred to a constrained 
field with a different mean. Near the roughening tem- 
perature, the switching procedure produces roughly 5% 
unswitched field points, and the corresponding mean field 
with these constraints can easily be determined using a 
steepest descent molecular dynamics method. To ensure 
ergodicity, a conventional Metropolis move is also carried 
out with every reverse MC move. 

Figure 1(a) shows simulation results for the scaling of 
the surface roughness a 2 with the length L of the lattice 
in simulations with different lattice sizes L 2 up to 1024 2 
and at several temperatures T from 16 to 30. KT theory 
26L 27 1 predicts a logarithmic divergence for a 2 with a 
universal slope at Tr 



a 2 (L) = *&(T R ) + -z]nL, 

7T 



(5) 



where a is the lattice constant of the surface, and in the 
units of Eqn. (4), a = 2ir. Therefore, at Tr the slope of 
Fig. 1(a) should be equal to 4. Above Tr, the logarith- 
mic behavior of a 2 continues to hold except the constant 
CTq as well as the slope both increase with T. The data 
in Fig. 1(a) show that for T — 21 and below, a 2 appears 
to approach a finite value as L — > oo. Therefore, it is 
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FIG. 2: Dynamic scaling for the relaxation time r (in units of 
MC passes) of a 2 as a function of lattice size L in Metropolis 
versus cluster MC, with their corresponding exponent £. 



clear that Tr > 21. The most recent simulation of the 
SG model by Sanchez et al. [32| (referred to as the "or- 
dered SG model" , OSGM, in this paper) suggested that 
Tr ss 16. Our data show that this is incorrect, and their 
error is likely due to slow sampling problems. Locating 
the precise value of Tr is more involved, since the data 
for T > 21 show no obvious tendency toward a finite a 2 . 
There are two possibilities: either these temperatures are 
above Tr or the system size may not be large enough to 
have reached the scaling limit for these temperatures. 
To determine which one is the case, we must resort to 
a comparison between the simulation data with KT the- 
ory. Figure 1(b) shows an expanded view of Fig. 1(a) 
for a few temperatures 23 < T < 30, but for each T 
the curve has been shifted vertically to remove the off- 
set ctq so that they all coincide at L = 64. The heavy 
dashed line indicates the KT slope at Tr according to 
Eqn.(5). The data therefore suggest that T R is slightly 
larger than 25 but less than 26, which is consistent with 
the RG prediction for Tr = 87r in the continuum model 
35L 3f|. The apparent lack of an asymptotic a 2 in the 
data for 21 < T < Tr implies that even for L = 1024, 
these lattice sizes are not yet large enough to be in the 
scaling limit for those temperatures. Finally, to compare 
the dynamic scaling behavior of the switching algorithm 
with Metropolis, Fig. 2 shows the relaxation time in the 
measurement of a 2 with the lattice size L slightly above 
Tr. Compared with the dynamic exponent £ w 2.5 in 
Metropolis, the switching algorithm shows a markedly 
improved £ ~ 1.4. 
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